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ABSTRACT 

The Fermi Gamma-ray Space Telescope reveals two large gamma-ray bubbles in the Galaxy, which 
extend about 50° (~ 10 kpc) above and below the Galactic center and are symmetric about the 
Galactic plane. Using axisymmetric hydrodynamic simulations with a self-consistent treatment of the 
dynamical cosmic ray (CR) - gas interaction, we show that the bubbles can be created with a recent 
AGN jet activity about 1-3 Myr ago, which was active for a duration of ~ 0.1 - 0.5 Myr. The bipolar 
jets were ejected into the Galactic halo along the rotation axis of the Galaxy. Near the Galactic 
center, the jets must be moderately light with a typical density contrast 0.001 < rj < 0.1 relative 
to the ambient hot gas. The jets are energetically dominated by kinetic energy, and over-pressured 
with either CR or thermal pressure which induces lateral jet expansion, creating fat CR bubbles as 
observed. The sharp edges of the bubbles imply that CR diffusion across the bubble surface is strongly 
suppressed. The jet activity induces a strong shock, which heats and compresses the ambient gas in 
the Galactic halo, potentially explaining the ROSAT X-idiy shell features surrounding the bubbles. 
The Fermi bubbles provide plausible evidence for a recent powerful AGN jet activity in our Galaxy, 
shedding new insights into the origin of the halo CR population and the channel through which massive 
black holes in disk galaxies release feedback energy during their growth. 

Subject headings: cosmic rays - galaxies: active - galaxies: jets - Galaxy: nucleus - gamma rays: 
galaxies 



1. INTRODUCTION 

Using 1.6 years of data from t he Fermi Gamma-ray 
Space Telescope^ Su et al.^ (l 2010l ) discovered two large 
gamma-ray bubbles emitting at 1 < < 100 GeV in 
our Galaxy. They are denoted as the 'Fermi bubbles', 
and extend to ~ 50° above and below the Galactic cen- 
ter (GC), with a width of about 40° in longitude. They 
have approximately uniform gamma ray surface bright- 
ness with sharp edges. The bubbles were previously iden- 
tified as the 'Fermi haze' from the first year Fermi data 
(jPobler et al.l[20Toh . The gamma-ray emission associ- 
ated with the bubbles has a hard spectrum {dN^/dE^ ~ 
E~'^) and could originate from cosmic ray (CR) pro- 
tons and/or CR electrons. CR protons may collide 
with thermal nuclei inelastically, producing neutral pi- 
ous w hich decay into gamma rays (jCrocker fc AharonianI 
l2011f ). For the latter, gamma rays can be produced 
by the inverse Compton (IC) scattering of CR electrons 
with the interstellar radiation field ( ISRF) and the cos- 
mic microwave background (CMB; ' Dobler et all I2010L 
iSu et al.l[20To[ ). The CR electron (CRe) population also 
produces microwave synchrotron radiation, which was 
previously invoked to explain the microwave haze (at 
tens of GHz) discovered by the Wilkinson Microwave 
Anisotropy Probe (WMAP) - the 'WMAP ha ze', which is 
spatially associated with the Fer mi bubbles (jFinkbeinerl 
[200l [PobTer fc Finkbeinerl [20081 ). 

The bilobular morphology of the Fermi bubbles sug- 
gests that the responsible CR population is not made 
of diffused CRs originally accelerated by supernova (SN) 
shocks in the Galactic plane. SN CRs diffusing into the 

^ UCO/Lick Observatory, Department of Astronomy and Astro- 
physics, University of California, Santa Cruz, CA 95064, USA; 
fulai@ucolick.org 



Galactic halo would likely have a much different spatial 
distribution, and tend to form one single CR-filled bub- 
ble centered at the GC. Furthermore, the sharp edges of 
the Fermi bubbles imply that the diffusion is strongly 
suppressed across the bubble surface, and thus is not the 
primary reason for the bubble expansion/formation (see 
Section [32]) • Similarly, the morphology and sharp edges 
of the Fermi bubbles suggest that the main CR popu- 
lation in the Fermi bubbles is not byproducts of dark 
matter annihilations. 

As pointed out by Su et al.l (|2010l ), the Fermi bub- 
bles were most likely created by some large episode of 
energy injection in the GC, e.g., jets originating from 
an active galactic nucleus (AGN) activity or a nuclear 
starburst. Here we perform the first dynamical study 
of the jet scenario. Our opposing jet model for the 
Fermi bubbles is inspired by similar non-thermal dou- 
ble jet events observed in massive galaxies, particu- 
larly tho se lyin^ at the centers of g alaxy groups and 
clusters (|McNamara fc NulsenI I2007D . AGN jets have 
been previously prop osed to explain extended extra- 
galactic radio sources (iLongair et an[197 3: Scheuer[ [l974l : 
iBlandfo rd fc Rees"1974y It is well known that AGN jets 
carry a significant amount of high energy CRs, which 
have been seen from radio synchrotron emission, and 
many AGN jets have also been observed to actively create 
CR-filled bubbles (e.g., the jets in Cygnus A and many 
Fanar off- Riley type I radio sources as in iLaing etld] 
2006). However, due to its limited sensitivity and resolu- 
tion, gamma-ray observation can not easily detect AGN 
bubbles. The proximity of the Fermi bubbles may pro- 
vide a special opportunity to thoroughly study AGN bub- 
bles in gamma ray bands. 

In this paper we use hydrodynamic simulations to per- 
form a feasibility study, investigating if the main mor- 
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phological features of the Fermi bubbles can be repro- 
duced by a recent AGN jet activity from the GC. Our 
calculations take into account the dynamical interactions 
between CRs and thermal gas, self-consistently model- 
ing the CR advection with the thermal gas. Our two- 
dimensional (2D) calculations allow us to do a fast pa- 
rameter study, identifying appropriate jet parameters. 
We show that to reproduce the morphology of Fermi 
bubbles, successful AGN jets cannot be either too mas- 
sive or too light, but are constrained to have densities 
about 0.001 - 0.1 times the density of the initial am- 
bient hot gas near the GC (i.e., with density contrasts 
0.001 <7^< 0.1). In particular, the observed Fermi bub- 
bles can be approximately reproduced by a pair of bipolar 
AGN jets which began about 1-3 Myr ago, and were 
active for a duration of 0.1 - 0.5 Myr. The total AGN 
energy released in our fiducial run Al is ~ 10^^ erg, but it 
scales with the uncertain density of the original gas in the 
Galactic halo, probably varying between about 10^^ and 
]^q57 YVe also show that to produce the sharp outer 
edges of the bubbles seen in gamma ray emission, CR 
diffusion across the bubble boundaries must be strongly 
suppressed below the CR diffusivity observed in the solar 
vicinity . 

Similar to other astronomical phenomena, the Fermi 
bubbles may result from the interplay of various physical 
mechanisms. Including hydrodynamics, CR dynamics, 
CR advection and diffusion, this paper mainly explores 
the potentially dominant physics involved in the forma- 
tion of the Fermi bubbles - a pair of opposing jets. We 
follow the dynamical evolution of the jets, and investi- 
gate the success and problems of the basic jet model by 
directly comparing with detailed observational features 
of the Fermi bubbles. We also investigate the degenera- 
cies and uncertaint ies associated w ith the jet model. In a 
companion paper ()Guo et ani2012t : hereafter denoted as 
"Paper 11" ) , we investigate the roles of additional physics 
- shear viscosity and CR diffusion on the evolution of the 
Fermi bubbles, exploring potential solutions for the ob- 
servational discrepancies associated with this basic jet 
scenario. The rest of the paper is organized as follows. 
In Section[2l we describe CR transport in the Galaxy, our 
Galactic model, and numerical setup. We present our re- 
sults in Section [3l and summarize our main conclusions 
with implications in Section IH 

2. THE MODEL AND NUMERICAL METHODS 

2.1. Cosmic Ray Transport in the Galaxy 

The observed Fermi bubbles are offset from the GC 
with one above and the other below the Galactic plane. 
They extend to about 50° 10 kpc) above and below 
the GC, and are symmetric about the Galactic plane. 
Here we assume that the distance from the Sun to the 
GC is Rq = 8.5 kpc. In Fermi maps, the bubbles are 
seen at 1 < < 100 GeV, corresponding to relativistic 
CRes with energies 10 < Ecr ^100 GeV if the observed 
7 rays are inverse Compton upscattering of starlight by 
these CRes. The IC cooling time of these CR electrons is 
typically a few million years, de pending; on thei r energies 
and locations (see Figure 28 of lSu et al.l[20T0h . Due to 
the short age and symmetry of these two bubbles, they 
probably share the same origin. The CRs may be ac- 
celerated or reaccelerated in situ in these two bubbles 



(jMertsch Sarkarll201lf ). but one potential problem for 
this scenario is to explain how turbulence is triggered 
roughly simultaneously in two bubble-shaped regions sep- 
arated by such a large distance (^ 10 kpc). Here we 
consider a different scenario where CRs are produced by 
one single event and then transported to these two bub- 
bles. If this scenario is correct, it is likely that the CRs 
were originated from the midpoint between the bubbles, 
the GC. 

Each CR particle has a velocity near the speed of light. 
However, it is well known that the collective transport 
speed of CRs in the Galaxy is much smaller (e.g., CRs 
reaching the Earth are observed to be highly isotropic). 
This is commonly explained by the scattering of CRs by 
magnetic irregularities. When the scattering is signif- 
icant, CRs are strongly trapped and 'effectively' move 
with magnetic irregularities, which are frozen in the hot 
plasma. Therefore, CRs are advected with thermal gas at 
the local gas velocity, and in this case, the CR transport 
is usually called advection. In regions where the mag- 
netic field is locally aligned, CRs may stream through 
thermal gas, but the streaming speed is limited by the 
Alfven speed. If CRs stream along the magnetic field 
lines at a speed faster than the local Alfven speed, the 
CR streaming instability is triggered and excite s Alfven 
wave s, which s ignificantly sc atter CRs in return (jSkillingI 
119711 : also see i KulsrudI 12005, Chap. 12). Consequently, 
the CR streaming speed relative to the gas is roughly 
the local Alfven speed. Due to the unknown and possi- 
bly complex magnetic field structure, CR streaming rel- 
ative to the local gas is often ignored, as we assume in 
this paper. This is a good approximation in many cases, 
where CR scattering off magnetic irregularities is signif- 
icant and/or the Alfven speed is much smaller than the 
local gas speed. 

In addition to advection, CRs can also diffuse through 
the thermal e; as, as they scatter off magnetic inhomo- 
geneities. CR diffusion is usually described by the CR 
diffusion coefficient a^, which probably depends on the 
magnetic field structure and CR energy. Typical values 
of f<i for 1 GeV CRs are found to h e k (3 - 5 ) x 10^^ 
cm^ s~^ in the solar vicinity (St rong et al.ll20Q7[ ). 

The Fermi bubbles extend up to ^ 10 kpc away from 
the GC. If the CRs in the bubbles are produced at 
the GC, they must be transported very quickly, at a 
speed of 'i;transport ~ 10 kpc/tage ^ 10000(tage/l Myr)"^ 
km/s, to form the bubbles, where tage is the age of 
the Fermi bubbles. We expect tage to be less than a 
few million years due to the short IC cooling times of 
CRes responsible for the detected 7-ray emission. Ob- 
viously, the main CR transport mechanism is not dif- 
fusion, which would form one single gamma-ray bub- 
ble centered at the GC and is too slow to transport 
CRs. To transport CRs for a distance of / = 10 kpc 
within tage by diffusion, the required diffusion coefficient 
is - ^ 3xl0^^(tage/l Myr)-^ cm^ s-\ three or- 

ders of magnitude larger than typical estimates of the CR 
diffusivity in the Galaxy. Furthermore, diffusion tends to 
produce blurred boundaries across which the CR density 
varies gradually, while the observed Fermi bubbles have 
very sharp edges in gamma rays. 

Since diffusion is not responsible for transporting CRs 
from the GC to the Fermi bubbles, advection of CRs 



Fermi BUBBLES CREATED WITH AGN JETS 



3 



with thermal gas is the only possible CR transport mech- 
anism. Strong CR advection may be induced by galactic 
winds or AGN jets from the GC. Galactic winds are often 
detected through absorption Hues of their cold gaseous 
components, which ty pically have speeds of about 200- 
300 km s~^ (Martin l2005|) , more than one order of 
magnitude below the required CR transport speed for 
the Fermi bubbles: 'i^transport ^ 10000(tage/l Myr)~-^ 
km/s. However, starburst winds usually contain multi- 
temperature components and some hot components may 
have much larger speeds (Strickland & Heckman 2009). 
Some winds have indeed been observed to have velocities 
slightly beyond 1000 km/s, particularly those probably 
drive n by radiation or mechan ical energy of AGN events 
(e.g., iRupke fc Veilleux|[2QTTI ). A thermally driven star- 
burst wind flowing at the speed of ^transport would re- 
quire gas with comparable sound speeds and tempera- 
tures T - 4 X lO^K that are unreasonably high. Alterna- 
tively, the wind scenario may be viable if the gamma- 
ray emission from the Fermi bubbles is mainly con- 
tributed by the decay of neutral pions produced during 
the hadronic collision s of CR pro tons with thermal nu- 
clei (Crocker & Ahar onianl l2011f ). One potential chal- 
lenge for the wind scenario is that winds often contain 
low- temperature gas as seen in Ha, molecular, and X-ray 
fllamentary structures, which have not yet been observed 
in the Fermi bubbles. 

In this paper we consider an alternative scenario where 
CRs are transported by AGN jets, which typically have 
relativistic or sub-relativistic velocities on parsec and kpc 
scales, much faster than galactic winds. The propagation 
of AGN jets is determined by the speed of hotspots at the 
jet extremities, and may be close to '^transport, depending 
on both jet parameters and the ambient gas properties. 
The primary goal of our paper is to study with numerical 
simulations the possibility of forming the Fermi bubbles 
with CR-carrying AGN jets from the GC. CRs in AGN 
jets may be dominated by the electron-proton plasma 
or electron-positron plasma. It is unclear if CR protons 
or electrons dominate the gamma-ray emission from the 
Fermi bubbles. In this paper we focus on the morphology 
and dynamical evolution of the Fermi bubbles, leaving 
studies of the particle content and gamma-ray emission 
mechanisms to the future. 

CR transport in the Galaxy has been studied 
numerically by some simulation codes, in partic- 
ular, the Galactic Propag ation (GALPROP) code 
(|Strong fc Moskalenkol Il998[ ). which is very detailed, 
including three-dimensional distributions of gas, dust, 
magnetic flelds, optical and FIR photons. However, 
GALPROP is not able to self-consistently model CR ad- 
vection, which is important when signiflcant gas motions 
are present. In addition to jet motions, CR pressure gra- 
dients can also produce gas motions, which in turn advect 
the CRs. Thus, dynamical interactions between CRs and 
the background gas should be taken into account to self- 
consistently model CR transport, particularly advection. 

We have developed a flnite-differencing code to numer- 
ically study CR transport and the dynamical interaction 
between CRs and hot gas, which has been successfully 
used to study X-ray cavities and radio bubbles in galaxy 
clusters in a series of papers, e .g., iMathews & Brighenti 
12008a'), 'Mathews & Brighenti' (2008b), Mathews (2009), 
IGuo Mathews. (2010a ). , Guo Mathews. (2010bl and 



IMathews fc Guol ()2010l ). More recently, we use our code 
to study the evol ution of CR-doni i nated AGN jets in 
galaxy clusters in i Guo fc Mathewsl ()2011f ) . In this pa- 
per, we intend to modify our code to study AGN jets 
and the formation of Fermi bubbles in the Milky Way. 
In the rest of this section, we will elaborate the basic 
equations and assumptions of our model, and the setup 
of our numerical procedure. 

2.2. Equations and Assumptions 

The dynamical interaction between CRs and thermal 
gas may be described by the CR pressure Pc. In our 
calculations, we directly follow the temporal and spatial 
evolution of the CR energy density ec, which is related 
to Pc through Pc = (7c — l)ec, where 7c = 4/3. The 
nature of the relativistic particles with energy density 
Cc is unspecifled and may be electrons and/or protons 
with any spectra. Of course the equation of state may 
be somewhat harder if Cc is mainly contributed by trans- 
relativistic protons at ~ 1 GeV. CR pressure gradients 
can accelerate thermal gas, converting CR energy into 
the kinetic energy of thermal gas, which may be further 
converted to the internal gas energy in shocks or other 
compressions. The combined hydrodynamic evolution of 
thermal gas and CRs can be described by the following 
four equations: 



1 + ^^-^ 



0, 



= -V(P + Pc) - pV$, 



de 

dt 



+ V ■ (ev) = -PV ■ V, 



dec „ , , 
^ + V.(e.v) 



-PcV-v + V-(/tVec), 



(1) 



(2) 



(3) 



(4) 



where d/dt = d/dt-\-v -S/ is the Lagrangian time deriva- 
tive, K is the CR diffusion coefficient, and all other 
variables have their usual meanings. The gas pressure 
P is related to the gas internal energy density e via 
P = (7 — l)e, where we assume 7 = 5/3. 

The Galactic potential <I> is mainly contributed by 
three components: the bulge, disk and dark matter halo. 
We assume that it is fixed with time and neglect the 
small contribution from hot halo gas. Our simple Galac- 
tic potential model and initial conditions for the thermal 
gas are explained in Section 12.31 We ignore radiative 
cooling of thermal gas, which is unimportant during the 
short-duration (< 1-3 Myr) of our simulations. The gas 
temperature is related to the gas pressure and density 
via the ideal gas law: 



T = 



kBP 



(5) 



where /cb is Boltzmann's constant, is the atomic mass 
unit, and /i = 0.61 is the molecular weight. To avoid 
confusion we denote the gas pressure P as Pg in the rest 
of the paper. 
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At time t = we assume that the CR energy density 
is zero in the Galaxy, ec = 0, but at later times CRs 
enter the Galactic halo in jets from the GC. Equation |4] 
describes the evolution of CR energy density including 
both advection and diffusion. CR diffusion depends on 
the magnetic field structure and may be anisotropic, but 
for simplicity, we generally assume that it is isotropic 
and the CR diffusion coefficient /i: is a constant in space 
and time: k ^ 3 x 10^^ - 3 x 10^^ cm^ s"^ (see Table 
[1]). We explore how the value of hi affects the formation 
and morphology of Fermi bubbles in Section 13. 2[ where 
we also consider a case with a spatially- varying diffusion 
coefficient. Hi may also depend on CR energy, which may 
affect the gamma-ray spectrum and the evolution of the 
Fermi bubbles. This effect is not studied in this paper. 

As discussed in the previous subsection, CRs interact 
with magnetic irregularities and Alfven waves, effectively 
exerting CR pressure gradients on the thermal gas (equa- 
tion [2]). Pressure gradients in the CR component act di- 
rectly on the gas by means of magnetic fields frozen into 
the gas. This is the primary interaction between CRs 
and thermal gas. We neglect other more complicated 
(probably secondary) interactions, e.g.. Coulomb inter- 
actions, hadronic collisions, and hydromagnetic- wave- 
mediated CR h eating, that all d epend on the CR energy 
spectrum (e.g. JGuo fc Oh|[2QQ8[ ). Since the jet evolution 
in our calculation is mainly driven by its kinetic energy 
and the injected CR energy, the dynamical effect of mag- 
netic fields is expected to be moderate and is unlikely to 
significantly alter our results (see our estimates in Sec- 
tion 13. ip . More sophisticated magnetohydrodynamical 
calculations are necessary in the future to explore the 
impact of magnetic fields. We also neglect CR energy 
losses from synchrotron and IC emissions. Although the 
cooling is important for TeV electrons, its impact on the 
integrated CR energy density Cc is expected to be small, 
i.e., the main contribution to Cc may be low-energy CR 
electrons (e.g., at ~ 0.1 - 100 GeV) and possibly CR 
protons, which have much longer lifetimes. If TeV elec- 
trons contribute significantly to the gamma-ray emission 
of the Fermi bubbles, CR cooling may have an impor- 
tant effect on the temporal evolution of the gamma ray 
emission (both intensity and spectrum), which is beyond 
the scope of the current paper. 



2.3. The Galactic Model 

Our key objective in this paper is to perform a fea- 
sibility study to determine if the Fermi bubbles can be 
created with jets from the GC that are apparently no 
longer active. The Fermi bubbles are inert ially confined 
by pre-existing gas in the Galactic halo having (virial) 
temperatures of ^ 10^ K. For simplicity, we assume that 
the hot gaseous halo is initially isothermal with tem- 
perature T = 2.4 X 10^ K, as inferred from X-ray ab- 
sorption line studies ()Yao fc Wangll2005f ). To derive the 
initial density distribution of the hot gas, we further as- 
sume that the gas is initially in hydrostati c equilibrium. 
We adopt the Galactic potential from iHelmi fc White! 
f2001), where the Galaxy is represented by a fixed po- 
tential with three components: a dark logarithmic halo 
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Fig. 1. — Initial profiles of electron number density of the hot gas 
in the Galactic halo along the R direction (parallel to the Galactic 
plane) at three fixed values of z (top panel) and along the z di- 
rection (perpendicular to the Galactic plane) at three fixed values 
of R (bottom panel) in all simulations except run A5, where the 
gas density is lower by a factor of 100. The density distribution 
ne(R, z) is solved from hydrostatic equilibrium, a spatially uniform 
temperature T = 2.4 x 10^ K, and a given value for the electron 
number density at the origin, rieo (chosen to be 1 cm~^ in all runs 
except run A5). 

a Miyamoto-Nagai disk 

^disc - 



GMdisc 



(7) 



and a spherical Hernquist stellar bulge 



bulge 



(8) 



where R is the galactocentric radius in the Galactic 
plane, r = \/ ^ z'^ is the spatial distance to the GC, 
^haio = 131.5 km s"^ and = 12 kpc; Mdisc = lO^^Mo, 
a = 6.5 kpc and 6 = 0.26 kpc; Mbuige = 3.4 x lO^^M© and 
= 0.7 kpc. We have also explored alternative models 
of t he Gala ctic potential from lBreitschwerdt et al.l ()1991l ) 
and Wolfire et al. (1995), and found that the resulting 
gas density distribution is similar. 

For the hot gas, it is convenient to use the electron 
number density ne, which is related to the gas density 
through 



(9) 



halo 



^halo 



ln(r^+d2), 



(6) 



where /ie = 5/i/(2 + /i) is the molecular weight per elec- 
tron. For a given value of electron number density at 
the origin neo, we derive the spatial distribution of Uq 
from hydrostatic equilibrium and the spatially-uniform 
temperature T = 2.4 x 10^ K. In all simulations except 
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run A5, we adopt the Uq distribution with neo = 1 cm~^, 
which gives electron number densities on the order of 
10~^ cm~^ at galactocentric distances of a few kpc. Fig- 
ure [T] shows the resulting Uq profiles along the R direction 
for three fixed values of z and along the z direction for 
three fixed values of R. We have also explored simula- 
tions with different values of neo, and found that, to re- 
produce the morphology of the observed Fermi bubbles, 
the required jet power scales with neo (i.e., with the den- 
sities of the hot halo gas; see discussions in Sec. 3.3). In 
run A5, we choose Uqo = 0.01 cm~^, explicitly exploring 
the dependence of our results on the density of halo gas. 
The unknown density distribution of hot halo gas is a 
major uncertainty in constraining the energetics of the 
Fermi bubble event. The observed gamma-ray surface 
brightness can in principle put a lower limit on the CR 
energy density within the bubbles, which may be used to 
constrain the properties of the hot halo gas and the local 
density of ISRF photons that are IC upscattered. 

2.4. Numerical Setup and Jet Injection 

Equations ([I]) — (jl]) are solved in (i?, z) cylindrical co- 
ordinates using a 2D axisymmetri c Eulerian code simi- 
lar to ZEUS 2D ( Sto nelT Norman! 1 1 992f ) . In particular, 
our code follows the evolution of both hot gas and CRs, 
implementing their dynamical interaction, CR diffusion, 
and an energy equation for CRs. The CR diffusion is 
solved by using operator-splitting and a fully implicit 
Crank-Nicolson differencing method. While this differ- 
encing scheme is stable, we restrict each time step by the 
stability condition for explicitly differenced diffusion, as 
well as Courant conditions for numerical stability. The 
computational grid consists of 400 equally spaced zones 
in both coordinates out to 20 kpc plus additional 100 
logarithmically-spaced zones out to 50 kpc. The central 
20 kpc in each direction is resolved to 0.05 kpc. For both 
thermal and CR fiuids, we adopt outfiow boundary con- 
ditions at the outer boundary and refiective boundary 
conditions at the inner boundary. 

The jet infiow is introduced with a constant speed '^jet 
along the z-axis (the rotation axis of the Galaxy) from 
the GC. During the active AGN phase < t < tjet, the 
jet is initialized in a cylinder on the computational grid 
with radius i^jet and length Zjet beginning at the GC. 
All the physical variables in this cylinder are updated 
uniformly with the jet parameters at every time step 
during this phase. The jet contains both thermal gas 
and relativistic CRs, which are initialized inside the jet 
source cy finder and injected into the Galaxy at z = zjet 
with the jet speed v^et having an initial opening angle 
of 0°. In the simulations presented in this paper, we 
take Zjet = 0.5 kpc. This allows us to avoid modeling 
the complex jet evolution closer to the GC, which would 
require higher spatial resolution and involves uncertain 
multi-phase gas properties. We further assume that the 
jet has undergone significant deceleration and transverse 
expansion within the jet source region (0 < z < Z]et), as 
seen in some extragalactic AGN jets (e.g., Laing et alj 
[2006). In our calculations, we thus investigate jets with 
RjQt/zjet ^ 0.2 — 1 and v^et/c ~ 0.05 — 0.3, where c is 
the speed of light. However, it is still possible that AGN 
jets from the GC may have much higher velocities at 
Zjet- Such jets are beyond the scope of this paper - our 
primary purpose is to perform a simple feasibility study 



to investigate if the location, size and morphology of the 
observed Fermi bubbles can be reproduced by a jet event 
within a few Myrs. 

We stop the calculation at the current age of the Fermi 
bubbles t = tpermi, which is defined as the time when the 
produced CR bubble reaches the most distant boundary 
of the observed northern Fermi bubble along the jet axis. 
The bubble age tpermi is constrained to be less than a few 
million years by the short IC cooling times of the gamma- 
ray-emitting CRes within the bubbles. In this paper we 
only consider models with tpermi ^ 3 Myr (the estimated 
IC cooling time of ^ 100 G eV electrons at z ~ 4 kpc; 
see Fig. 28 of Su et al.ll2010[ ), though it may be dynami- 
cally possible to create the Fermi bubbles with a current 
age more than 3 Myr. Alternatively, if the gamma-ray 
emission from the Fermi bubbles is mainly contributed 
by CR-proton-induced pion decay, the a^e of the Fermi 
bubbles could be much longer (see C rocker fc AharonianI 

ImiD. 

The initial jet is described by six parameters: the 
thermal gas density pj, the gas energy density ej, the 
CR energy density ejcr, the jet velocity Vjet^ the jet ra- 
dius i^jet, and the jet duration t^et- From pj, one can 
derive the commonly-used jet density contrast between 
thermal gas inside the jet and the ambient hot gas: 
V = Pj/Pamb, where Pamb IS the initial halo gas density 
at (i?, z) — (0, ^^jet). The jet power can be written as 

Pjet=Pke+Pcr+Pth, (10) 

where Pke = Pj'^fet^^j^et/^ the jet kinetic power. Per = 
ejcr'^jetTri^j^gt is the jet CR power, and Pth = eji'jetTri^j^g^ is 
the jet thermal power. The total energy injected by the 
jet can be written as E'jet = Pjet^jet- Thermal and CR 
energy density are related to the gas and CR pressure 
respectively, thus both affecting the dynamical evolution 
of hot gas. Therefore ej and ejcr are degenerate in our 
simulations: we can increase one while decreasing the 
other, reaching similar results (see run A6; however, syn- 
chrotron and/or IC emissions from CRs can in principle 
break this degeneracy) . In all simulation runs except run 
A6 (see Table[T]), we take r/e = ej/eamb = 1, where eamb is 
the initial ambient gas energy density at (i?, z) = (0, ^^jet). 
In simulations that successfully produce Fermi bubbles 
as observed, the jets are usually over-pressured relative 
to the ambient hot gas at z = Zjet- As explained above, 
the jet simulations have a very large parameter space, 
and our 2D calculations permit us to run a large number 
(more than 100) simulations with different parameters 
within a reasonable amount of time. Here we present 
a few representative calculations particularly relevant to 
the Fermi bubbles, and list their model parameters in 
Table [H Some basic information of the AGN event and 
the resulting Fermi bubbles in each simulation is given 
in Table El 

3. RESULTS 

3.1. Forming Fermi Bubbles with AGN Jets 

In this subsection, we present in detail one representa- 
tive run (run Al), which approximately reproduces the 
morphology of the Fermi bubbles as observed in gamma- 
rays with the Fermi telescope. The main objective is to 
show that it is feasible to produce the observed Fermi 
bubbles from one single AGN jet activity from the GC. 
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TABLE 1 

Model Parameters in All Simulations Presented in This Paper 









^jet 


^jet 


V 




^ d 
cjcr 


Vjet 


riei 


Run 


(cm ) 


(cm^ s""*") 


(kpcj 


(Myrj 










[LU cm ) 


Al 


1 


3 X 10^^ 


0.4 


0.3 


0.01 


1 


1.0 


O.lc 


5.68 


/\-QlIIl 


1 
1 


O X lU 


n A 


U.O 


U.Ul 


1 
1 


1 n 

l.U 


n 1 /- 


O.DO 


A-diff2 .... 


1 


3 X 10^9 


0.4 


0.3 


0.01 


1 


1.0 


O.lc 


5.68 


A-diff3 .... 


1 


varied 


0.4 


0.3 


0.01 


1 


1.0 


O.lc 


5.68 


A2 


1 


3 X 10^^ 


0.4 


0.3 


0.02 


1 


1.5 


O.lc 


11.36 


A3 


1 


3 X 1027 


0.4 


0.2 


0.01 


1 


3.0 


0.2c 


5.68 


A4 


1 


3 X 10^7 


0.2 


0.3 


0.05 


1 


6.0 


O.lc 


28.40 


A5 


0.01 


3 X 10^'^ 


0.4 


0.3 


0.01 


1 


0.01 


O.lc 


0.057 


A6 


1 


3 X 10^7 


0.4 


0.3 


0.01 


9.325 


0.1 


O.lc 


5.68 


Bl 


1 


3 X 1027 


0.4 


0.3 


0.0001 


1 


1.0 


O.lc 


0.057 


B2 


1 


3 X 1027 


0.4 


0.3 


0.5 


1 


1.0 


0.05c 


284 



^The initial thermal electron number density at the origin (R, z)= (0,0), which determines the 
density normalization of the isothermal hydrostatic halo gas (see Sec. 12. 3|) . 

^T] determines the initial thermal gas density in the jet base pj = Typatm, where Pamb is the initial 
halo gas density at (R, z) = (0, zjet) and zjet is chosen to be 0.5 kpc in this paper (see Sec. 12. 4|) . 

^?7e determines the initial thermal energy density in the jet base ej = ?7eeatm- 

^The initial CR energy density in the jet base (in units of 10 ^ erg cm ^). 

^riej is the initial thermal electron number density in the jet base: rigj = pj/(fiem^), calculated 
from the values of 77 and patm- 



TABLE 2 

Properties of the Fermi Bubbles and the Associated AGN Event in Simulation Runs 





^Fermi 


p a 
-t cr 








Mbu' 


AMbh^ 




Run 


(Myr) 






(1043 grg s-1) 


(1056 grg) 


(Mq/jt) 


(Mq) 


(Mo) 


Al 


2.06 


1.43 


7.09 


8.60 


8.13 


0.03 


9.0 X 10^ 


1.5 X 10^ 


A-diffl 


1.94 


1.43 


7.09 


8.60 


8.13 


0.03 


9.0 X 10^ 


1.5 X 10^ 


A-difr2 


1.30 


1.43 


7.09 


8.60 


8.13 


0.03 


9.0 X 10^ 


1.5 X 10^ 


A-difr3 


2.06 


1.43 


7.09 


8.60 


8.13 


0.03 


9.0 X 10^ 


1.5 X 10^ 


A2 


1.74 


2.15 


14.18 


16.41 


15.51 


0.057 


1.7 X 10^ 


3.0 X 10^ 


A3 


0.86 


8.58 


56.75 


65.48 


41.25 


0.23 


6.9 X 10^ 


2.0 X 10^ 


A4 


2.34 


2.15 


8.87 


11.03 


10.42 


0.039 


1.2 X 10^ 


1.9 X 10^ 


A5 


2.06 


0.014 


0.071 


0.086 


0.081 


0.0003 


90 


1.5 X 103 


A6 


2.15 


0.14 


7.09 


7.96 


7.52 


0.028 


8.4 X 10^ 


1.5 X 10^ 


Bl 




1.43 


0.07 


1.58 


1.49 


0.0055 


1.7 X 10^ 


1.5 X 10^ 


B2 


0.89 


0.72 


44.33 


45.08 


42.6 


0.16 


4.8 X 10^ 


3.7 X 10^ 



^Pcr, Pke, and Pjet are, respectively, the jet CR, kinetic, and total powers (in units of 10^^ erg s 
Pjet = Pke + Per + Pfch? whcrc the thermal jet power Pth is much smaller than P^e and/or Per in our runs. 
The real jet power depends on the uncertain halo gas density, and is thus not well constrained. 

^E'jet = Pjet^jet is the energy injected by one jet during the AGN phase < t < tjet- The total energy 
injected by both bipolar jets is 2^jet- 

^Mbh is the corresponding accretion rate of the supermassive black hole at the GC, assuming a feedback 
efficiency of 10%: Mrh = 2Pjet/(0.1c2). 

■^AMbh is the total mass accreted by the supermassive black hole at the GC during the AGN event, 
assuming a feedback efficiency of 10%: AMbh = ^BH^jet- 

^Mjet = 2pj7ri?2g^'Ujet^jet5 estimated at \z\ = 0.5 kpc, is the total mass ejected in the two jets during the 
AGN event, including the potential gas entrained by the jets at \z\ < 0.5 kpc. 



Other choices of jet parameters may produce similar re- 
sults, so it is impossible to determine unique jet prop- 
erties from the observed bubble morphology. Here we 
focus on some interesting features and problems in this 
representative run. 

The edges of the bubbles are shown in Figure 18 of 
iSu et al.l (j2QlQl ) in Galactic coordinates To com- 

pare our simulation results with observations, we derive 
the edges of the bubbles in our coordinates (i?, z) using 
the following coordinate conversion: 



cos{l) 



tan(6). 



(12) 



i? = i?0|tan(O|, 



(11) 



The resulting edge of the northern Fermi bubble is shown 
as the dotted region in each panel of Figure O The 
computed bubble is more vertically elongated in (i?, z) 
coordinates than when viewed in Galactic coordinates 
since dz/db ~ l/cos^(6) increases significantly with b 
from b = to b ^ 50°. The line-of-sight projection has a 
non-negligible effect on the coordinate conversion, since 
the bubble size is comparable to the distance between 
the Sun and the GC. The real bubble may be slightly 
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Fig. 2. — Central slices (16 x 15 kpc) of CR energy density (top panels) and thermal electron number density (bottom panels) in 
logarithmic scale in run Al at t = 1 Myr (left panels), and t = tpermi = 2.06 Myr (right panels). Horizontal and vertical axes refer to R and 
z respectively, labeled in kpc. The dotted region in each panel approximately encloses the observed north Fermi bubble. The propagation 
of the AGN jet, active for only tjet = 0.3 Myr, produces a CR bubble at t = 2.06 Myr approximately matching the observed Fermi bubble. 
The dashed lines in bottom panels trace the outer edge of the ROSAT X-ray shell feature in the northeastern direction (which is most 
prominent), and is roughly spatially coincident with the jet-induced shock at t = 2.06 Myr. 



smaller than the one shown in Figure 2. An accurate ac- 
count for this effect requires three-dimensional structure 
of the bubbles, and is much more complex. We expect 
this effect to be small and neglect it in this paper. 

The jet parameters for the fiducial run Al are listed 
below: radius i?jet = 0.4 kpc, thermal electron num- 
ber density Uq 



5.68 X 10- 



sity ej = 5.4 x 10 erg cm 



= 1.0 X 10 



gas energy den- 
CR energy density 
-9 cm~^, velocity Vjet = 0.1c, and 



cm 

-3 



duration tjet = 0.3 Myr. The jet is light with a density 
contrast = 0.01, which has been previously adopted in 
many jet simulations. Due to an additional CR compo- 
nent, the jet is over-pressured by a factor of ~ 10 with 
regard to the ambient hot gas at the jet base z = zjet- 
Energetically the kinetic power in the jet dominates over 
the CR power by a factor of ~ 5 (see Table [2j). The 
simulation is started at t = 0, and stopped at the cur- 
rent age of the Fermi bubbles t = tpermi = 2.06 Myr. 
The small value of tjet/^Fermi ~ 0.15 implies that all the 
CRs within the bubbles at the current time have nearly 
the same age, consistent with the Fermi observation that 
the gamma-ray spectr al index is quite uniform across the 
whole Fermi bubbles (jSu et al.ll20Tol ). 



Figure [2] shows central slices of CR energy density (top 
panels) and thermal electron number density (bottom 
panels) in logarithmic scale at t = 1 Myr (left panels), 
and t = tpermi = 2.06 Myr (right panels). As clearly 
seen, the jet produces a rapidly expanding CR bubble, 
which roughly matches the current north Fermi bubble 
(dotted circle) at t = 2.06 Myr. The CR transport is 
dominated by advection due to the high velocity of ther- 
mal gas in the jet. The CR diffusivity in this run is low 

= 3 X 10^^ cm^ s~^, having a very small effect, which 
can also be seen by the sharp edge of the resulting CR 
bubble (see Section [X2]) . The sharp edge is also an impor- 
tant feature of the observed Fermi bubbles. The lateral 
expansion of the bubble is mainly due to its high internal 
pressure (contributed by both CRs and shocked thermal 
gas) within the bubble and the rapidly decreasing am- 
bient gas pressure with the galactocentric distance. The 
resulting CR bubble is much narrower if the initial jet 
has a smaller internal pressure. Due to the expansion, 
the gas density in the bubble is significantly lower than 
the ambient gas density, as clearly seen in the bottom 
panels of Figure [2l This explains why ROSAT X-ray ob- 
servations show a 'cavity' of X-rays toward the center of 
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Fig. 3. — Variations of electron number density (top), the z- 
component gas velocity (middle), and pressures (bottom) along 
the jet axis for run Al at t = 1 Myr (dashed) and t = 2.06 Myr 
(dotted). The initial gas density and pressure profiles along the 
z-axis at t = are plotted as solid lines in the top and bottom 
panels, respectively. 
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Fig. 4. — Variations of electron number density (top), the R- 
component gas velocity (middle), and pressures (bottom) along the 
i?-direction (perpendicular to the jet axis) for run Al at t = 2.06 
Myr dX z = 2 (dotted), 5 (dashed), and 9 kpc (dot-dashed). The 
initial gas density and pressure profiles along the i?-direction at 
z = b kpc are plotted as solid lines in the top and bottom panels, 
respectively. 

the Fermi bubbles (similar to X-ray cavities often found 
in galaxy clusters). 

The bottom panels of Figure[2]show a strong shock pro- 
duced by the jet activity. The shock propagates into the 
hot Galactic halo, and at t = trermi (i-e., the time when 



the bubbles are currently being observed), has a speed 
of ~ 2000 km/s (mach number M ~ 9) at the north- 
most shock front. The propagating shock can also be 
seen in Figure [3l which shows variations of thermal elec- 
tron number density, the z-component of the gas velocity 
Vz, and pressures along the jet axis at t = 1 Myr and 
t = 2.06 Myr. In directions perpendicular to the jet axis, 
the shock has a slightly smaller speed, reaching R ^ 6 
kpc at t = tpermi, as sccu in Figure [21 Combining with 
the shock propagating toward the south direction, the 
Fermi bubble event produced a dumbbell-shaped shock 
front. The shocked gas near the shock front has higher 
densities {uq 1 - 4 x 10~^ cm~^) and temperatures 
(T ~ 3 - 7 keV), potentially explaining the dumbbell- 
shaped X-ray s hell features detected by the ROSAT 
X-ray telescope (|Sofuel [20001 : [Bland-Hawthorn fc CohenI 
12003). The ROSAT X-i ay sheh feature is most promi- 
nent in the no rtheastern d irection (the North Polar Spur; 
see Fig. 1 in [Sofue[[2000[ ), and its outer edge shown as 
the dashed line in Figure [2] is roughly spatially coincident 
with the shock front in the same direction. At t = tpermi, 
the shock front in the East direction has a speed of 

- 1800 km/s. At a distance of Rq/cos{35°) ^ 10.4 
kpc, the shock front (i.e., the outer edge of the X-ray 
feature) requires about 27 yrs to move an arcsecond, the 
approximate resolution of the Chandra X-ray Telescope. 
Our jet model can not reproduce the bi-conical struc- 
ture in the ROSAT 1.5 keV map as s hown in Figure 
5(d) of [Bland-Hawthorn Cohed (2003), which is lo- 
cated within the Fermi bubbles and may be produced 
by a more recent nuclear event. 

The strong shock (M ~ 8 - 9) shown in Figures 2, 
3 and 4 indicates that the Fermi bubbles are expanding 
explosively, only slightly decelerated by the small inertial 
resistance of the surrounding halo gas. The average total 
pressure inside our computed Fermi bubbles is about 50 

- 100 times larger than that in the un-shocked ambient 
halo gas. This is very different compared to radio bubbles 
and X-ray cavities in galaxy clusters, where usually weak 
shocks with Mach number M ~ 1 - 2 are induced. Fur- 
thermore, buoyancy which is widely recognized to play a 
key role in the evolution of X-ray cavities in many galaxy 
clusters, is much less important in the evolution of our 
simulated Fermi bubbles, where the momentum and en- 
ergy injected by the jets play the dominant role. 

The internal structure of the CR bubble can be seen 
in Figure [H which shows variations of electron number 
density, the i?-component gas velocity vr, and pressures 
along the i?-direction (perpendicular to the jet axis) at 
t = 2.06 Myr at three different heights: z = 2 (dotted), 
5 (dashed), and 9 kpc (dot-dashed). The bubble cor- 
responds to regions with high CR pressure, and is sur- 
rounded by the outward-propagating shock represented 
by large jumps in rig, ^'r and Pg, which show the struc- 
ture of the shock visible in Figure [2l The bubble has low 
gas densities ne ~ 1 - 8 x 10~^ cm~^, and is separated 
from the surrounding shocked gas through a contact dis- 
continuity across which the total pressure is continuous. 
Within the bubble, the CR pressure in run Al is typically 
~ 3 - 6 X 10~^^ dyn cm~^, which is comparable to the gas 
pressure, and slightly larger than the magneti c press ure 
of the 5 - 10 /iC fields considered bv [Su et al.[ (|2010[ ) to 
explain the "WMAP haze" observed at the base of the 
Fermi bubbles as synchrotron emission. This supports 
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our neglect of the magnetic field in the bubble dynam- 
ics, which will likely not significantly change our results. 
However, as discussed in the last paragraph of Section 
3.3, the CR pressure in our simulated bubbles depends 
on a model parameter - the initial CR pressure at the 
jet base, and does not represent the real CR pressure 
in the Fermi bubbles, which may be constrained by mi- 
crowave and gamma ray emissions from the bubbles (see 
Sec. 3.5). 

In run Al, the jet is energetically dominated by the 
kinetic energy, with the kinetic power Pke ~ 7.09 x 10^^ 
erg s~^ and the CR power Per ^ 1.43 x 10^^ erg s~^. 
The total jet power is Pjet ~ 8.6 X 10"^^ erg s'^. The 
jet is assumed to be produced by Sgr A*, the supermas- 
sive black hole located at the GC. Taking a standard 
accretion efficiency of 10%, the black hole accretion rate 
can be determined Mbh = 2Pjet/(0.1c^) = 0.03 Mq/jt 
(note that the black hole ejects two opposing jets). Dur- 
ing this AGN event {0 < t < t^et = 0.3 Myr), the black 
hole accreted a total of AMbh ^ 9000 Mq gas, and 
ejected in the jets a total of Mjet ~ 1.5 x 10^ Mq gas, 
which is estimated at \z\ = 0.5 kpc and thus includes 
the potential gas entrained by the jets at \z\ < 0.5 kpc 
(see Table 2). Assuming that Sgr A* has a mass of 
Mbh 4 X 10^ (iGhez et al.ll2008i) . the Eddington 
luminosity is PEdd 5.5 x 10"^^ erg s~-^, and the Edding- 
ton ratio for the AGN activity is e = 2Pjet/pEdd ~ 0.31. 
However, as we discuss in Section 3.3, the energetics of 
the Fermi bubble event scales with the uncertain density 
of the ambient halo gas, which confines the bubbles. For 
example, if the normalization of the halo gas is lower by 
a factor of 10, the energy and power of the associated 
AGN event are also smaller by the same factor 10. 

While the overall energetics, size, shape and age of our 
Al Fermi bubble shown in the right panels of Figure 2 
are encouraging, this bubble is deficient in two important 
respects that can be expected from any calculation based 
on the standard hydrodynamic equations 1-4. First, 
the pagoda-shaped bubble boundaries visible in Figure 2 
- formed by a series of Kelvin-Helmholtz vortices - dis- 
agrees significantly with the globally smooth, egg-shaped 
gamma-ray bub bles observed with the Fermi telescope 
(|Su et aLll2010f ). Second, the quasi-uniform distribution 
of CR electrons, proportional to ec{R^ z) plotted in Fig- 
ure 2, is inconsistent with the rather uniform gamma-ray 
surface brightness observed by Fermi across the bubble 
surfaces. The gamma-ray surface brightness produced 
by IC upscattering of rather smoothly distributed ISRF 
photons is roughly proportional to the integral of Cc along 
the line of sight through the bubble. For spatially uni- 
form ec, the surface brighness is expected to decrease 
significantly from the bubble center toward the bound- 
aries at constant z where the ISRF photon density is 
roughly constant or decreases slightly in the same di- 
rection. Since the ISRF photon density is expected to 
decrease in the z-direction (away from the Galactic cen- 
ter), ec{R^z) must also increase in this direction. Con- 
sequently, to produce a gamma-ray image with uniform 
surface brightness, it is necessary that ec{R^ z) increase 
smoothly toward all bubble boundaries. As we discuss 
below, these two observational discrepancies suggest that 
additional physics plays a significant role during the jet 
evolution. In Paper II, we show that these discrepancies 



can be corrected by including gas viscosity. 

3.2. Suppression of CR Diffusion across the Bubble 
Surface 

In our fiducial run Al, we choose a constant CR dif- 
fusion coefficient k = 3 x 10^^ cm^ s~^, which is much 
lower than typical estimates of diffusivity in our Galaxy 

^ (3— 5) X 10^^ cm^ s~^. A low value for hz is specifically 
chosen to produce the CR bubble with sharp edge, which 
is one of the main features in the observed Fermi bub- 
bles. When the diffusivity is not strongly suppressed, the 
edges of the resulting CR bubble are usually very broad, 
inconsistent with observations. To study the effect of CR 
diffusion on the bubble evolution, we perform three ad- 
ditional runs A-diffl, A-diff2, and A-diff3, which are all 
the same as run Al, but with different prescriptions for 
the CR diffusivity hz. 

In run A-diffl and A-diff2, the CR diffusion coefficient 
is chosen to be />: = 3 x 10^^, and 3 x lO^^cm^ s~^, re- 
spectively. In these cases, diffusion helps transport CRs, 
and thus the age of the CR bubble, tpermi = 1.94, 1.30 
Myr respectively, is shorter than that in run Al. The 2D 
distribution of CR energy density at t = tpermi in these 
runs is shown in left panels of Figure [5l When compared 
with run Al (right-top panel), CR diffusion in runs A- 
diffl and A-diff2 transports CRs to much large distances 
in the R direction, which is also clearly seen in Figure 
[6l which shows variations of CR pressure Pe along the 
P-direction at t = tpermi at z = 2 (bottom), and 9 kpc 
(top) in runs Al, A-diffl, A-diff2, and A-diff3. 

An important feature of the CR bubble produced in 
runs A-diffl and A-diff2 is that the bubble edge is not 
sharp, as clearly seen in Figures [5] and [6l This feature is 
clearly produced by CR diffusion, which transports CRs 
across the bubble surface, rendering the gradual outward 
decline in CR energy density at bubble edges. When the 
CR diffusion coefficient is significantly suppressed as in 
our fiducial run Al, the CR bubble is mainly expanded 
due to CR advection, and has sharp edges. Observa- 
tions of Fermi b ubbles show that the bubble edge is very 
sharp (Su et al. 2010), suggesting that the CR diffusion 
is significantly suppressed. However, the CR diffusion co- 
efficient may depend on the structure of local magnetic 
fields and may vary with space and time. One interesting 
question to ask is whether the sharpness of bubble edges 
implies that the CR diffusion is suppressed everywhere in 
the bubble regions (e.g., as in run Al). It is possible that 
the CR diffusion inside the bubble is not suppressed, as 
it does not directly affect the bubble surface. The key 
to produce sharp bubble edges is the suppression of CR 
diffusion across the bubble surface, which directly affects 
the structure of bubble edges. Physically, due to the 
small gyro-radii of relativistic particles (around 10~^ kpc 
for a 10 GeV electron in a 4 jaG field), CR diffusion may 
be anisotropic and CRs diffuse mainly along magnetic 
field lines with cross-field diffusion strongly suppressed. 
During the creation of the CR bubbles, ambient gas just 
external to the bubbles is strongly compressed, resulting 
in tangential magnetic field lines near the surface. This 
may explain why CR diffusion is suppressed across the 
bubble surface. It will be of great interest in the future 
to study this issue in detail using magnet ohydro dynamic 
simulations with anisotropic CR diffusion. 

Without explicitly introducing magnetic fields, here we 
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Fig. 5. — Central slices (16 x 15 kpc) of CR energy density in logarithmic scale in run A-diffl (top-left), Al (top-right), A-diff2 (bottom- 
left), and A-diff3 (bottom right) at t = ^Fermi (details listed in Tables [T] and [2] for each run). Horizontal and vertical axes refer to R and z 
respectively, labeled in kpc. The dotted region in each panel encloses the observed north Fermi bubble. In runs A-diffl and A-diff2, CR 
diffusion significantly affects the bubble evolution, rendering bubble edges that are less sharp than those observed. Variable CR diffusion 
in run A-diff3 leads to a smoother CR energy density distribution inside the bubble, while still suppressing CR diffusion across the bubble 
surface and retaining sharp bubble edges as in run Al. 

the bubble, since there are essentiahy no CRs there. At 
t = tpermi, the CR energy density distribution, shown 
in the bottom-right panel of Figure [5l is very similar 
to that in run Al (top-right panel), and particularly, 
the edges of the CR bubble are also very sharp, indi- 
cating that the prescription for CR diffusivity shown in 
equation [T3l indeed significantly suppresses CR diffusion 
across the bubble surface, and CR diffusion in the bub- 
ble interior is not required to be suppressed to produce 
the sharpness of the bubble edges. The bottom-right 
panel of Figure [5] also shows that the Cc distribution in- 
side the bubble in run A-diff3 is much smoother than in 
run Al. CR diffusion inside the bubble removes local 
CR structures (e.g., regions with high or low CR energy 
densities as seen in the right-top panel of Fig. 5), which 
may otherwise have been seen in the Fermi observations 
of projected gamma-ray emission. The observed Fermi 
bubbles have approximately uniform surface brightness, 
which may imply that CR diffusion is not strongly sup- 
pressed inside the bubbles (i.e., only the CR diffusion 
across the bubble surface is strongly suppressed). Fu- 
ture data from even longer-duration Fermi observations 
are needed to study the possible internal structure of the 
Fermi bubbles. 



want to do a preliminary study investigating how the 
CR bubble evolves if the CR diffusion is only suppressed 
across the bubble surface. To this end, we perform an 
additional run A-diff3, where the CR diffusion is normal 
{k, ~ 10^^ - 10^^ cm^ s~^) within the evolving CR bubble, 
but significantly suppressed exterior to the bubble sur- 
face. As shown in the CR bubble is separated from 
the surrounding thermal gas through a contact disconti- 
nuity, across which thermal gas density changes abruptly. 
Thus we assume that in run A-diff3, the CR diffusivity is 
related to thermal gas density through an ad-hoc equa- 
tion: 

^ _ fa X 10^^(neo/ne) cm^ s~^ when Ue > ^eo^^) 
\3 X 10^^ cm^ s~^ when Uq < neo, 

where Uqq = 10~^ cm~^. The parameters in equation 
[T3l are chosen so that during the calculation of this run 
(t < tpermi), CR diffusion is always significantly sup- 
pressed outside the expanding bubble {k, < 10^^ cm^ 
s~-^), but not suppressed within it {hz ~ 10^^ - 10^^ cm^ 
s~^). The low CR diffusivity outside the bubble only 
suppresses CR diffusion across the bubble surface and 
does not directly affect regions much further away from 
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Fig. 6. — Variations of CR pressure Pc along the i?-direction at 
t = tFermi ^^,t z = 2 (bottom), and 9 kpc (top) in runs Al (solid), 
A-diffl (dotted), A-diff2 (dashed), and A-difTS (dot-dashed). The 
edge of the CR bubble is not sharp in runs A-diffl and A-diff2 
because of CR diffusion across the bubble boundary. 

3.3. A Parameter Study and Degeneracies 

The jet simulation involves a large number of parame- 
ters, many of which are degenerate. It is thus impossible 
to uniquely determine all jet properties from the mor- 
phology of the observed Fermi bubbles. For example, a 
shorter jet duration tjet results in a longer bubble age and 
a 'fatter' bubble, but this effect can be averted by a larger 
jet density nej, a higher jet velocity 'Ujet, or a smaller CR 
energy density ejcr in the jet. Here we present results of 
a few additional (essentially degenerate) runs (A2, A3, 
and A4), all of which roughly produce the morphology of 
Fermi bubbles. The jet parameters and some results are 
summarized in Tables [T] and [21 Compared to run Al, the 
jet is more massive in run A2 (77 = 0.02), faster in run A3 
('^jet = 0.2c), and narrower in run A4 {R^et = 0.2 kpc). 
The age of Fermi bubbles in these runs is tFermi = 1-74, 
0.86, 2.34 Myr respectively. The 2D distribution of CR 
energy density at t = tpermi in these runs is shown in 
Figure [71 where the electron number density distribution 
in run A2 is also shown in the bottom-left panel. In 
run A2, the shock fits even better with the outer edge 
of the X-ray shell feature (dashed line) than in run Al, 
strengthening the point that the ROSAT X-idiy features 
surrounding the Fermi bubbles are produced by shocked 
gas associated with the jet activity. These four runs are 
representative of successful jet models which roughly re- 
produce the overall size and shape of the Fermi bubbles, 
suggesting that the Fermi bubbles could be produced by 
a powerful AGN jet event which began about 1-3 Myr 
ago, and was active for a duration of ~ 0.1 - 0.5 Myr. 

The evolution of AGN jets and the resulting CR bub- 
bles also depend on the initial density of the confining 
halo gas. For a gaseous halo with fixed temperature, its 
initial hydrostatic density distribution can be scaled up 
or down by varying the central gas density neo- Equa- 
tions [H - [H ensure that the evolution of the jet and re- 
sulting CR bubble is the same if the following three jet 
properties - density pj, thermal energy density ej and 



CR energy density ejcr are scaled by the same factor as 
UeQ. As an example, here we present the result of an ad- 
ditional run A5, where neo, Pj, ej, and ejcr are all scaled 
down by the same factor 100 with respect to run Al. 
The spatial distribution of CR energy density in this run 
at time t = tpermi is shown in the top panel of Figure 
[HI confirming that the morphological evolution of the jet 
and bubble is the same as in run Al. It is important 
to note that the jet duration and the current age of the 
Fermi bubbles are not affected by the halo gas density, 
while the total jet energy and the CR energy density 
within the bubbles scale with it. Thus, due to our limited 
knowledge of the halo gas density, the energetics of the 
AGN event responsible for the Fermi bubbles - as well 
as the CR energy density and the gamma-ray emissivity 

- is not constrained very well by our model. Depending 
on the initial density distribution in the Galactic halo, 
the energetics of the Fermi bubble event may be ~ 10^^ 

- 10^'^ erg (Table 2). 

In principle, accurate X-ray observations of the shell 
of shocked gas surrounding the Fermi bubbles can im- 
prove our estimate of the gas density distribution in the 
pre-bubble atmosphere and consequently the energetics 
of the dynamical event that created the bubbles. All- 
sky RO SATX - idiy ob servations at 0.75 and 1.5 keV were 
used bv iSofuel (|2000l ) to identify dumbbell-shaped shells 
of soft X-ray emission having lobes about 10 kpc in di- 
ameter symmetric about the Galactic center and aligned 
approximately along the Galactic rotation axis. At the 
time Sofue (2000) interpreted this as emission from gas 
in the Galactic halo that was shocked by an explosion 
in the Galactic center with total energy £^tot ^ 10^^ erg 
that occurred about 1.5 x 1 0^ yr a^o. H owever the orig- 
inal atmosphere adopted by ISofuel ((2000) is considerably 
different (not in hydrostatic equililDrium) and about 18 
times less dense than ours at = 0, z = 5 kpc (Figure [T]) 
near the bubble c enter. The co rresponding single-bubble 
energy found by ISofuel (|2000l ). when scaled to our at- 
mosphere, is Ejet = 0.5^tot X 18 ^ 9 X 10^^ erg, which 
is consistent with our model Al where Ejet = 8 x 10^^ 
erg. Similar to the Fermi bubbles, the real energetics of 
the ROSAT X-ray shells is not well constrained due to 
the uncertain density of the original gas in the Galactic 
halo. However, if they were indeed produced by a jet 
event (e.g. the Fermi bubble event), their age may be 
much shorter than 15 Myrs as estimated in a starburst 
model bv Sofue (2000). 

The evolution of AGN jets and resulting bubbles de- 
pends on the total pressure in the initial jet, but it is 
insensitive to how the jet pressure is distributed between 
thermal and non- thermal components. Thus, another 
uncertainty in our dynamical model is the ratio of the 
CR to thermal pressure within the jet, which is directly 
related to an important problem - the particle content 
(thermal versus non-thermal) within the current Fermi 
bubbles. In run Al, the CR pressure dominates in the 
initial AGN jet (the ratio of CR to thermal pressure is 
9.25), resulting in comparable CR and thermal pressures 
within the Fermi bubbles at the current time (see the 
bottom panel in Fig. [4]). However, the CR-to-thermal 
pressure ratio in the bubbles depends on its initial value 
at the jet base, and is generally not 1 in our simulations. 
In an additional run A6, we follow the evolution of a jet 
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Central slices (16 x 15 kpc) of log(ec) in run A2 (top-left), A3 (top-right), A4 (bottom-right), and log(ne) in run A2 (bottom 
(details listed in Tables [T] and [5] for each run). Horizontal and vertical axes refer to R and z respectively, labeled in kpc. 



Fig. 7. 

left) at t = tpei 

The dotted region in each panel encloses the observed north Fermi bubble. The dashed lines in the bottom-left panel trace the outer edge 
of the ROSAT ^-ray emission feature surrounding the northern bubble. Run Al and these additional runs have different jet parameters, 
but all produce CR bubbles approximately matching the observed Fermi bubbles. 



dominated by thermal pressure (the initial ratio of CR 
to thermal pressure is ~ 0.1; the total jet pressure, jet 
density and other parameters are the same as run Al; see 
Table [T]) , finding that the morphological evolution of the 
resulting CR bubbles in this run is very similar to run Al, 
as clearly seen in the bottom panel of Figure [51 Due to 
much lower CR power in the initial jet, run A6 produces 
Fermi bubbles having proportionally lower CR pressure 
Pc (much lower than thermal pressure) and lower gamma 
ray surface brightness than our fiducial run Al. The CR 
pressure in our successful runs by no means represents 
the real CR pressure in the Fermi bubbles, which can 
be constrained by microwave and gamma ray emissions 
from the Fermi bubbles (see Sec. 3.5). However, our dy- 
namical model is robust and the bubble evolution, which 
depends on the total jet pressure, is insensitive to the 
CR pressure. Analyses based on multi- wavelength (X- 
ray, gamma-ray, microwave, radio) observations of the 
Fermi bubbles and the surrounding gas, which is beyond 
the scope of this paper, should be used to constrain the 
density of halo gas and the particle content of the Fermi 
bubbles. 

3.4. Constraints on the Jet Density 



As discussed in the previous subsection, it is hard to 
accurately determine the jet properties, but it is possible 
to put some useful constraints, particularly on the jet 
density. When the jet is very light, e.g., r] = 0.0001 as in 
run Bl (see Table [T] for other jet parameters), it quickly 
decelerates in the hot gas, and expands laterally in the 
R direction due to the high pressure within the bubble. 
The top panel of Figure [9] shows the CR energy density 
in logarithmic scale at t = 3 Myr. As clearly seen, the jet 
only reaches z ^ 6 kpc at this time, and the resulting CR 
bubble is quasi-spherical, unlike the observed Fermi bub- 
bles elongated in the z direction. Very light jets tend to 
have strong backfiows (i.e., with Vz < inside the bubble 
at large i?), which expand laterally, for ming; CR bubbles 
much 'fatter' than Fermi bubbles (see iGuo fc Mathewsl 
l2011f ) . We have experimented with more runs with differ- 
ent jet densities, and find that the Fermi bubbles can be 
successfully reproduced only with jet densities rj > 0.001. 

On the other hand, when the jet is very massive, e.g., 
T] = 0.5 as in run B2 (see the bottom panel of Figure [9]), 
it decelerates very slowly, producing very little backfiow. 
In this case, most of the CRs are advected by the jet to 
the jet tip, not forming a 'fat', radially-elongated bubble 
as observed. We have run many simulations with high 
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run A5, t = 2.06 Myr 




run B1 , t = 3 Myr 
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run A6, t = 2.15 Myr 



Fig. 8. — Central slices (16x 15 kpc) of log(ec) in run A5 (top) and 
A6 (bottom) at t = tFermi- Horizontal and vertical axes refer to R 
and z respectively, labeled in kpc. The dotted region in each panel 
encloses the observed north Fermi bubble. Run A5 shows that the 
required jet power and the CR energy density are proportionally 
smaller in a gaseous halo with lower densities compared to run 
Al, while run A6 shows that the morphological bubble evolution 
is quite insensitive to the relative contribution of CRs and thermal 
gas to the jet pressure when the total jet pressure and other model 
parameters are fixed. 

jet densities, and found that the formed CR bubbles are 
usually much thinner than the observed Fermi bubbles. 
Acceptable runs typically have the jet density contrast 

??<o.i. 

AGN jets efficiently transport CRs from the GC to the 
Galactic halo. Additionally the jets can produce or reac- 
celerate CRs in the strong shock at the jet tip (hot spot) 
which we have not considered here. The morphology 
and evolution of the resulting CR bubble significantly de- 
pend on many jet properties, but the jet density contrast 
is particularly important. From the morphology of the 
Fermi bubbles, we can approximately constrain the jet 
density contrast: 0.001 ^ ^ ^ 0.1, which corresponds to 
6 X 10~^ cm~^ ^ ^ej ^ 6 X 10~^ cm~^ for our adopted 
model for the Galactic thermal atmosphere. The con- 
straint on the jet density contrast r] is quite robust, not 
depending on the value of neo, the density normalization 
of the ambient halo gas. 

3.5. Success and Problems of Our Jet Model 

The detection of the Fermi bubbles in the Galaxy is 
one of the biggest discoveries in astronomy in recent 
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Fig. 9. — Central slices (16 x 15 kpc) of CR energy density in 
logarithmic scale in run Bl at t = 3 Myr (top) and B2 (bot- 
tom) at t = ^Fermi- Horizontal and vertical axes refer to R and 
z respectively, labeled in kpc. The dotted region in each panel en- 
closes the observed north Fermi bubble. Very light jets in run 
Bl (?7 = 0.0001) decelerate rapidly in the halo gas, forming quasi- 
spherical CR bubbles unlike the radially elongated Fermi bubbles 
observed. On the other hand, the massive jets in run B2 (77 = 0.5) 
transport most CRs just along the z axis, not forming the 'fat' 
bubbles observed. 

years. The origin of the bubbles remains mysterious. 
Here we are proposing that the bubbles were produced 
by a recent powerful AGN jet event from the GC. We 
focus on the formation and dynamical evolution of the 
Fermi bubbles, leaving details of the multi- wavelength 
emission features to future work. In this subsection, we 
summarize the success and potential problems of our jet 
model in explaining the observed features of the Fermi 
bubbles. 

Our model successfully explains a few key observa- 
tional features of the Fermi bubbles. (1) The origin of 
CRs is naturally explained within a jet scenario. CRs can 
be accelerated to very high energies in AGN jets near su- 
permassive black holes and/or in jet hotspots, as seen in 
numerous observations of extragalactic AGN jets. (2) 
Sub-relativistic jets naturally produce large CR bubbles 
within a few Myrs, a short bubble age inferred from the 
fact that the gamma-ray emission is li kely dominate d by 
the I C emission of CRes (jPobler et al.ll2Q10,. and Su et al.l 
I2OIOI ). It is not trivial to transport or accelerate CRs 
in 10 kpc sized bubbles within such a short amount of 
time. (2) The opposing jet scenario naturally explains 
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the bilobular morphology of the bubbles, which are sym- 
metric with respect to the GC. (3) We identified a few 
sets of jet parameters, with which the jets produce CR 
bubbles with morphologies similar to the observed bub- 
bles (mainly the bubble elongation and ellipticity). (5) 
In our model, the jet duration is much shorter than the 
bubble age, implying that the CRs in the whole bubble 
have similar age at the current time. This explains the 
Fermi observation that the gamma ray spectral index is 
quite uniform across most regions of the two bubbles. 
(6) The jet event produces a strong shock, which heats 
and compresses the ambient gas in the Galactic halo, 
potentially explaining the ROSAT X-ray shell features 
surrounding the bubbles. (7) Sharp bubble edges indi- 
cate that the Fermi bubbles have been expanding mainly 
due to CR advection, while CR diffusion across bubble 
edges is strongly suppressed. 

However, there are still some potential problems asso- 
ciated with our model, which need to be investigated in 
future studies. The main difference between the simu- 
lated CR bubble and the observed Fermi bubbles is that 
the former suffers from Kelvin-Helmholtz instabilities at 
its surface (as clearly seen in Figure [2j), while the real 
Fermi bubbles have smooth edges. Line-of-sight projec- 
tion tends to smooth the bubble edge when viewed in 
Galactic coordinates, but it is unable to completely re- 
move surface irregularities created by these instabilities, 
as seen in Figure [TOl which shows for run Al the line- 
of-sight projection of CR energy density (a proxy for the 
gamma-ray surface brightness due to CR electrons) and 
pec (a proxy for the gamma-ray surface brightness origi- 
nated from CR protons). This discrepancy suggests that 
there may be additional physics, which plays an impor- 
tant role in the bubble evolution, suppressing the devel- 
opment of these instabilities. Similar surface instabilities 
have received much attention in numerical studies of ra- 
dio bubbles a nd X-ray cavities in ga laxy clusters, where 
gas v iscosity (iRevnolds et al.l [2QQ5[ ) and magnetic ten- 
sion (Ruszkowski et ani2QQ7[ ). have been invoked to suc- 
cessfully suppress these instabilities at the boundaries of 
buoyantly-rising bubbles. But it remains to be shown 
that the instabilities can also be suppressed during the 
creation of AGN bubbles. A preliminary study exploring 
the beneficial role of viscosity on the bubble evolution 
and morphology is presented in Paper II. Viscous sup- 
pression of instabilities in the Fermi bubbles may pro- 
vide an unusual opportunity to study the interesting gas 
microphysics - viscosity - in hot plasma. 

An interesting feature of the observed Fermi bubbles 
is the approximately uniform gamma-ray surface bright- 
ness, particularly at high latitudes {\b\ > 30°). How- 
ever, a nearly uniform CR distribution, as in our mod- 
els computed here, produces a center-brightened (limb- 
darkened) surface brightness distribution in projection 
(see Figure [TQ|) . which is inconsistent with uniform 
gamma ray brightness observed in the Fermi bubbles. 
A fiat surface brightness distribution implies that the 
CR density increases gradually toward bubble edges and 
with Galactic latitude from the bubble center. However, 
if CRs are too concentrated at the bubble surface, the 
gamma-ray surface brightness will be too strongly limb- 
brightened. Thus, the CR distribution may need some 
fine-tuning to get a roughly fiat gamma-ray intensity. 
Although such an edge-favored CR distribution is not 
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Fig. 10. — Line-of-sight projection of CR energy density (top; a 
proxy for the gamma-ray surface brightness due to CR electrons) 
and pcc (bottom; a proxy for the gamma-ray surface brightness 
originated from CR protons) in logarithmic scale in run Al at 
t = 2.06 Myr. Horizontal and vertical axes refer to Galactic longi- 
tude and latitude respectively, labeled in degrees. The solid circle 
represents the edge of the observed south Fermi bubble. Edge ir- 
regularities and edge darkening are clearly seen at high latitudes. 
At low latitudes (|6| < 20° — 30°), Fermi observations are sig- 
nificantly contaminated by emissions from the Galactic disk and 
bulge. The spatial variations of CR spectra and ISRF may affect 
the gamma-ray map, which is not considered here. 

reproduced in our current jet simulations, this discrep- 
ancy does not invalidate the jet scenario for the Fermi 
bubbles. Instead, it may require a new physical mecha- 
nism that operates additionally in the standard AGN jet 
model. In Paper II, we argue that hot gas viscosity may 
provide the additional physics to suppress surface irregu- 
larities and to produce an edge-favored CR distribution, 
resulting in a more uniform gamma ray surface brightness 
distribution. Other physical mechanisms (e.g., stochastic 
re-acceleration of CR electrons preferentially ne ar bubble 
edges as suggested bv lMertsch fc Sarkarll2QlTI ) may also 
contribute to surface brightness uniformity. 

In this paper, we choose a very simple model for the 
ambient gas distribution - a hot volume-filling isother- 
mal gas in hydrostatic equilibrium. In reality, the gas 
distribution in the Galaxy may be much more complex. 
Gas motions (e.g. winds) may be present, producing non- 
axisymetric features in the bubble morphology. At low 
latitude, the hot, warm and cold gas within and near the 
Galactic bulge may significantly affect the jet evolution 
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Fig. 11. — Gamma-ray spectrum from a CR electron popula- 
tion (IC scattering; solid line) and a CR proton population (pion 
decay; dashed line). The CRe spectrum is taken to be a power 
law in energy: dN/dE oc with 0.1 < E < 1000 GeV. The 

CR proton spectrum is assumed to be a power-law in momentum: 
dN/dp oc p~'^-^ with 0.1 < p/rripC < 1000. Crosses represe nt the 
obser ved average gamma ray flux from the Fermi bubbles ( Su et al.l 
l20T0h . 

and the bubble shape there. Some of our current simu- 
lations produce a 'stem' feature at low latitude (see Fig. 
[2j), which can not be accurately studied without model- 
ing the multi-phase gas there. Due to strong gamma-ray 
contaminations from the Galactic disk and bulge at low 
latitude, further studies are required to investigate if the 
stem is indeed present or not in the observed Fermi bub- 
bles. 

With hydrodynamic simulations, we focus on the evo- 
lution and morphology of the Fermi bubbles in this pa- 
per. Our simulations also investigate the transport of 
CRs (advection and diffusion) by following the evolution 
of the CR energy density ec, which may be contributed by 
CR electrons and/or protons with arbitrary spectra. The 
real CR content and spectrum in the Fermi bubbles may 
be constrained by the observed gamma ray spectrum, 
which is be in^ actively studied and debated in the liter - 
ature (e.g. ISu et all [20101 : ICrocker fc AharonianI IMTI ) . 
Here we present simple emission calculations in two ex- 
treme cases, assuming that the gamma ray flux is emit- 
ted entirely by CR electrons and protons r espectively. In 
the leptonic scenario, we follow ISu et al.l (^010), adopt- 
ing a power-law electron spectrum dN/dE oc E~'^-^ with 
0.1 < £^ < 1000 GeV. We choose the normalization of 
the electron spectrum so that the calculated gamma ray 
spectrum, which is shown as the solid line in Figure [TTJ 
is roughly consistent with the observed gamma ray lumi- 
nosity. Here we used the K lein-Nishina IC cross section 
ffilume nthal fc Gouldll 19701 ) and an approximate line-of- 
sight integration length of 4 kpc across the bubbles. The 
ISRF model is taken from the latest version (version 
54.1.984) of GALPROP (evaluated at the bubble cen- 
ter = 0, |z| = 5 kpc). We also calculated synchrotron 
emission from CR electrons in this model and found that, 
when the magnetic field strength in the bubbles is cho- 
sen to be 8 /iC, the synchrotron flux at 23 GHz is around 
1.1 kJy/sr, consistent with the detected WMAP flux at 
around h ~ —20° - —30°. It remains unclear why the 
microwave flux fro m the south bubble drops significantly 
at 6 < -35° (see Dobled 120121 ). 

Combining with our simulations, emission calculations 
may provide new insights on the Fermi bubbles. To this 
end, we calculated the CR pressure contributed by elec- 
trons between 0.1 and 1000 GeV in the above leptonic 



model and found that it is around = 2.0 x 10~^^ dyn 
cm~^, which is about four orders of magnitude below the 
total pressure (~ 10~^^ dyn cm~^) inside the simulated 
Fermi bubbles in run Al (see Figs [3] and 3]). This result is 
quite robust with respect to the power index of the elec- 
tron spectrum, suggesting that the pressure within the 
Fermi bubbles is not dominated by CR electrons in the 
leptonic scenario. The dominant bubble pressure may 
instead come from thermal gas (e.g. in run A6), CR 
protons, or magnetic fields. 

In the hadronic scenario, the gamma ray emission of 
the Fermi bubbles is dominated by the decay of neutral 
pions produced from hadronic collisions of CR p rotons 
with thermal nuclei ([Crocker fc AharonianI |2011[ ) . We 
assume that CR protons have a power-law distribution 
in momentum: dN/dp cx with 0.1 < p/rrinC < 1000. 
We use an analytic formula in'EnBlin et ajj (|2007l eq. 62) 
to approximate gamma-ray emissivity from pion-decay 
and adopt thermal electron density in the Fermi bub- 
bles to be 10""^ cm~^ (as in run Al; Figs [3] and |4]), 
which was used to derive the target nucleon density. The 
resulting gamma-ray flux is shown in Figure [11] (dot- 
ted line) and the CR proton pressure in this model is 
Pc = 5.1 X 10~^^ dyn cm~^, slightly larger than the to- 
tal pressure (~ 10~^^ dyn cm~^) inside the simulated 
Fermi bubbles in run Al. This does not rule out the 
hadronic scenario, since (1) the required CR proton pres- 
sure drops if the thermal gas density in the bubbles is 
higher and (2) the total bubble pressure in our dynamical 
simulations increases if the initial ambient gas pressure 
is larger. However, our simple calculations suggest that 
the required CR proton pressure in the hadronic scenario 
may dominate the total pressure in the bubbles, and is 
much larger than the CR electron pressure in the leptonic 
scenario. 

4. SUMMARY AND IMPLICATION 

The detection of Fermi bubbles in the Galaxy is one 
of the most important and striking findings during the 
first two years' observations of the Fermi Gamma-ray 
Space Telescope. The bubbles are symmetric about the 
GC, with one above and the other below the Galactic 
plane. The surface brightness in 7-ray emission from the 
bubbles is quite uniform with sharp edges at the bub- 
ble boundaries. These observed properties make it very 
difficult to explain the bubbles in many of the standard 
ways, for example by supernova shocks in the Galactic 
plane or by dark matter annihilations. 

The unique location, morphology, and sharp gamma- 
ray edges of the bubbles suggest that they were proba- 
bly created by an episode of energy injection in the GC. 
Gamma rays from the bubbles can be produced by CR 
electrons and/or protons. Motivated by numerous dou- 
ble radio sources and X-ray cavities observed in other 
galaxies, we investigate in this paper if AGN jets from 
the GC can create the Fermi bubbles with the observed 
morphology, size and appropriate age. In our model, a 
pair of bipolar jets were released from the GC along the 
rotation axis of the Galaxy. We model the jet evolution 
using a series of 2D axisymmetric hydrodynamical simu- 
lations, following the evolution of CRs and their dynam- 
ical interactions with thermal gas (see Section [2?2]) . CR 
advection is modeled self-consistently and CR diflFusion 
is also considered. 
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We show that the observed Fermi bubbles could be 
reproduced by a recent AGN jet event about 1-3 Myr 
ago, which was active for a duration of ~ 0.1 - 0.5 Myr. 
The total jet energy released during the AGN event may 
be ~ 10^^ - 10^^ erg, depending on the initial density dis- 
tribution of the hot halo gas which confines the bubbles. 
Containing both thermal gas and CRs, the two-fluid jet 
quickly advects CRs to high latitude. The jet is energeti- 
cally dominated by the kinetic energy, and over-pressured 
with either CR or thermal pressure which induces lateral 
expansion, creating a fat CR bubble as observed. In our 
fiducial run Al, the two opposing jets have a total power 
of 2Pjet ~ 1.7 X 10^^ erg s~^, corresponding to an ac- 
cretion rate of Mbh ~ 0.03 Mq/jt for the central black 
hole and an Eddington ratio of e ~ 0.31. The jet activity 
also induces a strong dumbbell-shaped shock (currently 
Mach number M ~ 8 - 9), which heats and compresses 
the ambient gas in the Galactic halo, qualitatively simi- 
lar to the dumbbell-shaped X-ray shell features observed 
by the ROSAT X-idij telescope. 

The two observed Fermi bubbles are roughly located 
and elongated along the rotation axis of the Galaxy. 
In the jet scenario, this morphology requires that the 
two opposing jets that produced the Fermi bubbles were 
ejected nearly along the Galactic rotation axis, as as- 
sumed in our simulations. This jet direction may be 
possible if the gas accreted by Sgr A* had angular mo- 
mentum in this direction or Sgr A* had a relatively 
strong spin along this direction before the accretion. 
It is unclear if this special orientation observed in the 
Fermi bubbles is a coincidence or a general feature of 
CR bubbles in disk galaxies (see Baum e tld] 11993 ' for 
the same preferential orientation detected in radio bub- 
bles in Seyfert galaxies). 

Because of the degeneracy of successful bubble models 
with different jet parameters, it is difficult or impossi- 
ble to determine unique jet parameters directly from the 
bubble location and morphology. However, we can put 
some useful constraints on them, particularly on the jet 
density contrast relative to the ambient hot gas at the jet 
base. Very low density, ultra-light jets decelerate quickly 
and usually form quasi-spherical bubbles unlike the ob- 
served Fermi bubbles which are radially elongated. In 
contrast, heavy, massive jets decelerate very slowly, ad- 
vecting most CRs to the tip of the jet with too little 
lateral expansion. Successful jets are moderately light, 
having typical density contrasts 0.001 ^ ^ ^ 0.1 (i.e., 
6 X 10~^ cm~^ ^ ^ej ^ 6 X 10~^ cm~^ for our adopted 
model for the Galactic thermal atmosphere). 

We also show that to produce sharp edges of the bub- 
bles, CR diffusion across the bubble surface must be 
suppressed significantly below the CR diffusion rate ob- 
served in the solar vicinity. CR advection is responsible 
for the spatial expansion of the Fermi bubbles, which 
compresses thermal gas near the bubble surface, aligning 
the local magnetic fields to be tangential to the bub- 
ble boundary. CR diffusion is probably anisotropic, with 
cross-field diffusion strongly suppressed. However, it is 
likely that CR diffusion deeper inside the bubbles is not 
suppressed and it tends to remove small CR structures, 
resulting in a smoother CR distribution and 7-ray emis- 
sion within the bubbles as suggested by Fermi observa- 
tions. 



The total pressure inside the Fermi Bubbles must not 
be less than the pressure of CRs required to produce the 
gamma ray emission observed. If the gamma-ray emis- 
sion from the bubbles is mainly due to CR electrons, the 
required CR electron pressure to produce the observed 
gamma-ray flux is negligible compared to the total bub- 
ble pressure, which may instead be dominated by other 
components, e.g., thermal gas, CR protons, or magnetic 
fields. On the other hand, if the gamma-ray emission 
is mainly due to CR protons, the required CR proton 
pressure is much higher, probably dominating the total 
bubble pressure. 

The Fermi bubbles provide plausible evidence for a 
recent powerful AGN jet activity in the Milky Way. 
Extragalactic jets have been detected in n iany distant 
galaxies , such as the famous jet in M87 0unor et al.l 
1999; Kovalev et al. 2007). Numerous bipolar radio bub- 
bles and X-ray cavities have also been observed, clearly 
associated with the central stellar bulges of massive 
galaxies. These non- thermal regions are almost cer- 
tainly produced by AGN jets, which are often not de- 
tectable because of the short time when they are active 
(McNamara & Nulsen 2007). Multi- wavelength studies 
show that AGN jet activity deposits large amounts of 
mechanical and CR energy into the host environments, 
significantly affecting the evolution of gas and host galax- 
ies. But in our Galaxy there has been little evidence for a 
currently- active AGN jet. It would be very unusual if Sgr 
A* has always been quiescent in the past, as during its 
growth, a large amount of energy O.IMbhc^ ^ 10^^ 
erg) may have been released. The detection of Fermi 
bubbles as the remnant of a recent powerful AGN jet 
event is thus of dramatic astronomical importance, sug- 
gesting that AGN jet activity may also happen regularly 
in our Galaxy. 

Furthermore, AGN activity from the GC may con- 
tribute significantly to, or dominate, the CR population 
in the Galactic halo. Suppose for example that new pairs 
of Fermi bubbles with total CR energy 2.7 x 10^^ ergs 
(twice that of the single bubble in run Al) are produced 
every 30 Myr (AGN duty cycle ^ 1%). This would sup- 
ply CR energy to the Galactic halo at a rate 2.86 x 10^^ 
ergs s~^. By comparison, if one supernova of energy 
IQ^i ergs occurs in the Galaxy every 50 years and 10% 
of that energy is converted to CRs, the total CR power 
generated, 0.6 x 10^^ ergs s~^, is considerably less than 
the estimated rate from AGN activity. Nevertheless, CRs 
from previous generations of Fermi bubbles may not con- 
tribute significantly to the locally observed CR energy 
density - this would require that CRs return to the so- 
lar vicinity by relatively slow diffusion in the halo before 
buoyancy and kinetic energy of the moving bubbles carry 
them entirely out of the Galactic halo. 

Do similar Fermi bubbles exist in other galaxies? Due 
to the limited sensitivity and resolution of gamma-ray 
observations, the Fermi telescope can not easily detect 
extragalactic Fermi bubbles similar to those in the Milky 
Way. However synchrotron emission of many bubbles 
has been detected and spatially resolved in radio obser- 
vations. AGN jets have been observed in many spiral 
Seyfert galaxies where extended k pc-scale radio s truc- 
tures (KSRs) are common ( Galli ml^7rd~d] [2006h . In 
a sample of 12 Seyfert KSRs, iBaum et all (jl993[ ) found 
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that the extended radio emission tends to ahgn with the 
projected minor axis of the host galaxy, similar to the ori- 
entation of t he Fermi bubble s with respect to the Galac- 
tic disk fsee lElmouttie et alJ [T998 and Kha rbet al.ll2QQ6l 
for two well-studied examples). It remains to be deter- 
mined if the Fermi bubbles and extragalactic KSRs are 
indeed similar astrophysical phenomena. 

Finally, as in Sections 3.1 and 3.5, we draw attention 
to two shortcomings of our CR+hydrodynamic computa- 
tions of the Fermi bubble evolution: surface irregularities 
and non-uniform gamma ray limb darkening. Irregular- 
ities in the bubble boundaries are induced by instabili- 
ties in the differentially shearing flow between vertically 
downflowing gas in the bubbles and the (assumed) ini- 
tially stationary gas in the Galactic halo. In addition, the 
rather uniform distribution of cosmic ray energy density 



inside our computed bubbles would produce IC gamma 
ray intensities that decrease toward the bubble bound- 
aries, assuming that the ISRF photons (which are up- 
scattered) are smoothly distributed. These inadequacies 
may suggest that a new physical mechanism in addition 
to the basic jet hydrodynamics, such as shear viscosity 
discussed in the companion paper - Paper II, plays a role 
during the evolution of the Fermi bubbles 
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